function varargout = prbfit (x, pp)

   ## usage:  varargout = prbfit (x, pp)
   ##
   ## estimate probit using bspline with k knots, extrap. with rate q

   global PAR

   if isempty(PAR)
      PAR.ncpu = 1 ;
      PAR.par_opt = {"UniformOutput", false, "VerboseLevel", 0} ;
   endif

   P = 1e-10 ; # strong spline smoothing

   [N, x] = twodim(x) ;
   [nr nc] = size(x) ;

   if nargin > 1
      CJ = {} ; for i=1:nc, CJ = {CJ{:}, {x(:,i), pp(:,i)}} ; endfor 
      w = parcellfun_(@(x) _prbfit(x), CJ, PAR.par_opt{:}) ;
      varargout{1} = reshape(cell2mat(w), N) ;
   else
      CJ = mat2cell(x, nr, ones(1, nc)) ;
      w = parcellfun_(@(x) _prbfit({x}), CJ, PAR.par_opt{:}) ;
      w = cell2mat(w) ;
      varargout{1} = reshape(w(1,:), 1, N(2:end)) ;
      varargout{2} = reshape(w(2,:), 1, N(2:end)) ;
   endif 

endfunction


function varargout = _prbfit (x)

   ## usage:  varargout = _prbfit (x)
   ##
   ## wrapper for prbfit

   global PAR PRSV ;

   ppo = PAR.ppo ;

   if Lin = (length(x) > 1), pp = x{2} ; endif 
   ## transition to new pp.form element
   if Lin && !isfield(pp, "form")
      pp.form = "pp" ;
      pp.breaks = pp.x(:)' ;
      pp.coefs = pp.P ;
      pp.pieces = size(pp.coefs, 2) ;
      pp.order = size(pp.coefs, 3) ;
      # pp.orient = 1 ;
      pp.dim = pp.d ;
   endif
   ## transition to new pp.form element
   x = x{1} ;

   P = 1e-10 ; # strong spline smoothing

   z = nan(rows(x), 1) ; 

   k = 5 ;

   if !any(I = !isnan(x)) && !Lin
       warning("xds:xds", "prbfit: no values") ;
       pp.breaks = nan(1,2*k-1) ;
       pp.dim = nan ;
       pp.coefs = nan(length(pp.breaks)-1,ppo) ;
       # pp.orient = nan ;
       pp.pieces = size(pp.coefs, 2) ;
       pp.order = size(pp.coefs, 3) ;
       pp.form = "pp" ;
       varargout{1}(1:2,1) = pp ;
       return ;
   endif

   if Lin
      if any(isnan(pp.breaks))
	 z(I) = nan ;
      else
	 z(I) = ppval(pp, x(I)) ;
      endif 
      varargout{1} = z ;
      return ;
   endif

   u = x(I) ;

   [F, IF] = cdf(u) ;
   if PRSV
      [wu um us] = anom(u) ;
      F = norminv(F, um, us) ;
      if PAR.dbg, printf("prbfit: %f %f\n", um, us) ; endif
   else
      F = norminv(F) ;
   endif
   if PAR.prbraw
      varargout = {nan, F} ;
      return ;
   endif 

   ## extrapolate and fit
   u = u(IF) ; F = F(IF) ;

   ub = getbreaks(u, k) ;
   Fb = getbreaks(F, k) ;

   uu = extrap(u) ;
   FF = extrap(F) ;

   pp = spfit(uu, FF, ub, ppo) ;
   ppinv = spfit(FF, uu, Fb, ppo) ;

   if PAR.dbg
      pp_print(u, F, pp, ppinv) ;
   endif

   ## transition to new pp.form element
   if !isfield(pp, "form")
      pp.form = "pp" ;
      pp.breaks = pp.x ;
      pp.coefs = pp.P ;
      pp.pieces = size(pp.P, 2) ;
      pp.order = size(pp.P, 3) ;
      pp.dim = pp.d ;
   endif
   if !isfield(ppinv, "form")
      ppinv.form = "pp" ;
      ppinv.breaks = ppinv.x ;
      ppinv.coefs = ppinv.P ;
      ppinv.pieces = size(ppinv.P, 2) ;
      ppinv.order = size(ppinv.P, 3) ;
      ppinv.dim = ppinv.d ;
   endif
   ## transition to new pp.form element

   varargout{1}(1,1) = pp ;
   varargout{1}(2,1) = ppinv ;
#   varargout{1} = pp = csaps(uu, FF, P) ;

   if (nargout > 1)
      z(I) = ppval(pp, x(I)) ;
      varargout{2} = z ;
   endif

endfunction


function xb = getbreaks (x, k, d)

   ## usage: xb = getbreaks (x, k, d)
   ## 
   ## get breakpoints

   if nargin < 3, d = 3 ; endif

   xm = (x(1) + x(end)) / 2 ;
   xd = (x(end) - x(1)) / 2 ;

   xd = linspace(0, xd^d, k).^(1/d) ;

   xb(k+1:2*k) = xm + xd ;
   xb(1:k) = xm - fliplr(xd) ;

endfunction


function xx = extrap (x, q)

   ## usage: xx = extrap (x, q)
   ## 
   ## extrapolate x on both sides

   if nargin < 2, q = 0.3 ; endif

   m = rows(x) ; kq = min(m-1, floor(m*q)) ;
   xx = nan(m + 2*kq, 1) ;

   xx(1+kq:m+kq) = x ;
   xx(1:kq) = 2*x(1) - flipud(x(1+1:1+kq)) ;
   xx(m+kq+1:m+kq+kq) = 2*x(m) - flipud(x(m-kq:m-1)) ;
	 
endfunction


function output = spfit(x,y,xb,ord)
%SPFIT Fit a spline to noisy data.
%   PP = SPFIT(X,Y,XB) fits a piecewise cubic spline with break
%   points XB to the noisy data (X,Y). Use PPVAL to evaluate PP.
%
%   Example:
%       x = 2*pi*(0:100)/100;
%       y = sin(x) + 0.1*randn(size(x));
%       xb = 0:6;
%       pp = spfit(x,y,xb);
%       z = ppval(pp,x);
%       plot(x,y,'.',x,z);
%
%   PP = SPFIT(X,Y,XB,N) is a generalization to piecewise polynomial 
%   functions of order N (degree N-1), with continuous derivatives
%   up to order N-2. Default is a cubic spline with N = 4.
%
%       N	Function
%       -----------------
%       1	Step function
%       2   Piecewise linear and continuous
%       3   Piecewise quadratic with continuous first derivative
%       4   Piecewice cubic with continuous second derivative (default)
%       5   Etc.

%   Author: jonas.lundgren@saabgroup.com, 2007.

if nargin < 1, help spfit, return, endif
if nargin < 2, y = 1; endif
if nargin < 3, xb = 0; endif
if nargin < 4, ord = 4; endif

%   Check data vectors

x = x(:);
y = y(:);
m = length(x);
if length(y) ~= m
    if length(y) == 1
        y = y*ones(size(x));
    else
        error('Data vectors x and y must have the same length!')
    endif
endif

%   Sort and check the break points

xb = sort(xb(:));
xb = xb([diff(xb)>0; true]);
n = length(xb) - 1;

if n < 1
    n = 1;
    xb = [min(x); max(x)];
    if xb(1) == xb(2)
        xb(2) = xb(1) + 1;
    endif
endif

%   Adjust limits

xlim = xb;
xlim(1) = -Inf;
xlim(end) = Inf;

%   Generate power- and coefficient-matrices for smoothness conditions

as = [ones(1,ord); ones(ord-1,1)*(ord-1:-1:0)-(0:ord-2)'*ones(1,ord)];
as = max(as,0);
cs = cumprod(as(1:ord-1,:));
ps = as(2:ord,:);
B0 = cs.*0.^ps;

%   Smoothness conditions

B = zeros((ord-1)*(n-1),ord*n);
h = diff(xb);
for k = 1:n-1
    Bk = cs.*h(k).^ps;
    B((ord-1)*(k-1)+1:(ord-1)*k, ord*(k-1)+1:ord*(k+1)) = [Bk, -B0];
endfor

%   QR-factorization

nn = min(size(B));
[Q,R] = qr(B');
Q2 = Q(:,nn+1:end);

%   Weak conditions (least square sense)

A = zeros(m,n+ord-1);
a = zeros(m,1);
mm = 0;
for k = 1:n
    I = (x <= xlim(k+1)) & (x > xlim(k));
    xdata = x(I) - xb(k);
    ydata = y(I);
    d = length(xdata);
    Ak = (xdata*ones(1,ord)).^(ones(d,1)*(ord-1:-1:0));
    A(mm+1:mm+d,:) = Ak*Q2(ord*(k-1)+1:ord*k,:);
    a(mm+1:mm+d) = ydata;
    mm = mm + d;
endfor

%   Solve

c = Q2*(A\a);

%   Make piecewise polynomial

coefs = reshape(c,ord,n).';
output = mkpp(xb,coefs);

return

endfunction


function pp_print (u, F, pp, ppinv, pfile)

   ## usage:  pp_print (u, F, pp, ppinv, pfile)
   ##
   ## print pp diagnostics

   persistent nprint = 0 ; nprint++ ;
   global PAR PRJ BATCH ;

   q = 2 ;

   if (nargin < 5) 
      pfile = fullfile("data", PRJ, PAR.pdd, ["var_" num2str(nprint) ".png"]) ;
   endif
   if BATCH
      graphics_toolkit gnuplot ; figure("visible", "off") ;
   endif

   ue = affin(u, q) ;
   Fe = affin(F, q) ;

   subplot(1,2,1) ; plot(u, F, ".", ue, feval(@ppval, pp, ue)) ;
   xlabel("phys.") ; ylabel("norm.") ;
   subplot(1,2,2) ; plot(F, u, ".", Fe, feval(@ppval, ppinv, Fe)) ;
   xlabel("norm.") ; ylabel("phys.") ;

   if BATCH
      print(pfile) ;
      close ;
      printf("pp_print: --> %s\n", pfile) ;
   endif

endfunction
